#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "eno.h"
#include "fluids.h"

#include "eno.xtn"
#include "fluids.xtn"


#if defined(CRAY)
/* the cray convention is to capitalize fortran names, including common
   block names
*/
#        define FF_xaverage XAVERAGE
#        define FF_yaverage YAVERAGE
#        define FF_pwrspec PWRSPEC
#else
/* the sun convention is to give fortran names a trailing underscore */
#        define FF_xaverage xaverage_
#        define FF_yaverage yaverage_
#        define FF_pwrspec pwrspec_
#endif
 
/*
   This program sucks up ascii files as created by amrplot and
   performs the indicated diagnostics
*/
#define INFILE_LEN 85

static float *data;
static float *xavg,*yavg;
static float *freq,*power,*specdata;
static float *xmarc,*ymarc;
static int numpoints,nx,ny;
static float xlen,ylen,dx,dy,cutoff;

static double wave_array[128];
FILE *wave_file;

double power_max;

void init_wavedata()
{
int i;

  wave_file=fopen("wavedata","r");
  for (i=0;i<128;i++)
   {
    fscanf(wave_file,"%lf\n",&(wave_array[i]) );
   }
  i=fclose(wave_file);
}
 
double get_wavedata(xval)
double xval;
{
double hw,x1,x2,f1,f2,f3;
int ix,ix1,ix2;

  hw=10.0/128.0;
  f1=(xval/hw)-0.5;
  ix=(int) (floor(f1));
  if (xval<hw/2.0) ix=-1;
  ix1=ix;if (ix1<0) ix1=127;if (ix1>127) ix1-=128;
  ix2=ix+1;if (ix2>127) ix2-=128;
  x1=(ix+0.5)*hw;x2=(ix+1.5)*hw;
  f1=wave_array[ix1];
  f2=wave_array[ix2];
  f3=f1+((f2-f1)/(x2-x1))*(xval-x1);
  return f3;
}

void init_marcus()
{
int i,j;
  
  nx=mx;
  ny=my;
  numpoints=nx*ny;
  xlen = 10.0;
  ylen = 5.0;
  dy = ylen/ny;
  dx = xlen/nx;
  cutoff=dx*2.0;


  data = (float *) calloc(numpoints,sizeof(float));
  xavg = (float *) calloc(ny,sizeof(float));
  yavg = (float *) calloc(nx,sizeof(float));
  ymarc = (float *) calloc(ny,sizeof(float));
  xmarc = (float *) calloc(nx,sizeof(float));
  /*  psd of y-averaged data  */
  freq = (float *) calloc(nx/2,sizeof(float));
  power = (float *) calloc(nx/2,sizeof(float));
  specdata = (float *) calloc(nx,sizeof(float));

  for(j=0; j<ny; j++) {
   *(ymarc + j) = dy*(j+0.5) - 0.5*ylen;
  }

  for(i=0; i<nx; i++) {
   *(xmarc + i) = dx*(i+0.5);
  }

}  /* init_marcus */

void my_pwrspec()
/* e.g. nx=64  dx=10/64  specdata:0..nx-1   power:0..nx-1 */
{
int nx_half,k,j;
double domain,sum_re,sum_im,beta;

  nx_half=nx/2;
  domain=dx*nx;
  for (k=0;k<nx_half;k++)
   {
    sum_re=sum_im=0.0;
    for (j=0;j<nx;j++)
     {
      beta=2.0*pi*k/domain;
      sum_re+=specdata[j]*cos(beta*j*dx);
      sum_im+=specdata[j]*sin(beta*j*dx);
     }
    power[k]=(sqr(sum_re)+sqr(sum_im));  /* no sqrt */
   }
}
    

void dump_marcus(time_step_count)
int time_step_count;
{
FILE *output1, *output2, *output3;
int i, j,ix,iy;
char outfile1[INFILE_LEN];
char outfile2[INFILE_LEN];
char outfile3[INFILE_LEN];

  sprintf(outfile3,"%s.psd_den",file_name);
      
  if (time_step_count==0) 
   output3 = fopen(outfile3,"w");
  else
   output3 = fopen(outfile3,"a");

/*    read data  */
  i=0;
  for_all_points 
   if (point_ok)
    {
     *(data+i)=gdt(1,ix,iy);
     i++;
    }
 
/*  x-average as a function of y  */
  FF_xaverage(&nx,&ny,&cutoff, data, xavg);
/*  y-average as a function of x  */
  FF_yaverage(&nx,&ny,&cutoff, &dy, data, yavg);

/*  psd of y-averaged data  */

  for(i=0; i<nx; i++) {
   *(specdata + i) = *(yavg + i);
  }
/*
  FF_pwrspec(&nx,&dx,specdata,freq,power);
*/
  my_pwrspec();
 
/*  dump  */
  if (time_step_count==0)
   {
    fprintf(output3,"%5d\n",(nx/2)-1); 

    power_max=0.0;
    for(i=1;i<nx/2;i++)
     power_max=max(power_max,fabs( *(power+i) ) );
    if (power_max<almost_same)
     power_max=1.0;
   }

  for(i = 1; i < nx/2; i++) {
   fprintf(output3," %17.10e\n", (*(power + i))/power_max); 
  }
printf("in dump_marcus\n");
  i=fclose(output3);
 
}
 
